RMarkdown history: 1. Feb 2021 started github repo (https://github.com/shu251/midcayman-rise-microeuk) 2. continue to commit changes and as you include knitr and commit changes, those will be pushed to the git repo 3. eventually set up the gitpages under settings 4 https://resources.github.com/whitepapers/github-and-rstudio/

1 Mid-Cayman Rise grazing

Summary of cruise operations

Below is the compilation of all microscopy data

1.1 Set up working environment

GitHub repo hosts all raw data to reproduce analysis and results. See input-data or clone the repo.

library(tidyverse); library(cowplot); library(broom)
getwd()
## [1] "/Users/sarahhu/Desktop/Projects/Mid-Cayman_Rise/midcayman-rise-microeuk"

2 Protistan cell concentration

Import eukaryotic cell count data from grazing experiments. In this section, we will calculate cells per ml from raw counts (Field of view, etc.) and use to estimate protist cell concentration. These will be used below in grazing experiment calculations.

counts <- read.delim("input-data/euk-counts-compiled.txt", 
                     blank.lines.skip = FALSE,
                     na.strings = c("", "NA"),
                     stringsAsFactors = FALSE) # Import
counts[is.na(counts)] <- 0 # Change blanks to zeroes

Raw data table collected during microscopy count process. Below code reviews the structure of this raw data and updates column headers to be more ‘R’ friendly.

colnames(counts) <-c("DATE", "SAMPLE", "EXPID", "VOL", "MAG", "FOV", "nanoNoFLP", "microNoFLP", "nanoFLP", "microFLP", "NOTES", "DateCompiled"); colnames(counts)
##  [1] "DATE"         "SAMPLE"       "EXPID"        "VOL"          "MAG"         
##  [6] "FOV"          "nanoNoFLP"    "microNoFLP"   "nanoFLP"      "microFLP"    
## [11] "NOTES"        "DateCompiled"
# Column 9 and 10 list the total number of FLP inside each, so # of occurences

# Notes if count was not countable
# unique(counts$NOTES)

To count occurence and number of FLP ingested by eukaryotic cells, the number of FLPs ingested was tallied and comma separated for multiple eukaryotic cells with FLP. These values need to separated and counted as 1 eukaryotic cell each, but retain the number of FLP per cell.

counts_occur <- counts %>%
  filter(NOTES != "Not countable") %>% 
  mutate(nanoFLP_occur = as.numeric(str_count(nanoFLP, "[1-9]\\d*")), #Count number of euk cells observed with FLPs
         microFLP_occur = as.numeric(str_count(microFLP, "[1-9]\\d*")),
         nanoTOTAL = as.numeric(nanoNoFLP) + nanoFLP_occur, #Add number of euk cells with FLPs to those without for total number of euk cells
         microTOTAL = as.numeric(microNoFLP) + microFLP_occur,
         euksTOTAL = nanoTOTAL + microTOTAL) %>%
      data.frame

2.1 Microscopy cells per ml (euk)

Input data are the raw microscopy counts by FOV. Code below calculations cells/ml based on these values. Additionaly, variance and standard deviation are also calculated. Eukaryotic cells were also classified by size, where micro equates to >20um and nano is <20um. All counts were done at 100x magnification, confirm this:

unique(counts_occur$MAG)
## [1] 100
# head(counts_occur)

Calculate cell concentration (cells/ml).

counts_cellsml_all <- counts_occur %>%
  group_by(SAMPLE, EXPID, VOL) %>% #Calculate averages by sample
  summarise(totalFOV = n(), # Count total FOV counted
            nanoAvg = sum(nanoTOTAL)/totalFOV, #Average per FOV
            nanoVar = var(nanoTOTAL), #Variance
            nanoSd = (2*(sqrt(nanoVar))), #Standard deviation
            microAvg = sum(microTOTAL)/totalFOV, ## Repeat for microeuks
            microVar = var(microTOTAL), 
            microSd = (2*(sqrt(microVar))), 
            euksAvg = sum(euksTOTAL)/totalFOV, ## Repeat for total cell count
            euksVar = var(euksTOTAL), 
            euksSd = (2*(sqrt(euksVar))), 
            .groups = 'drop_last') %>%
  # Calculate cells/ml based on magnification (at x100, 0.01 is vol of grid), volume filtered (VOL), dilution factor (0.9), and area of counting grid (for Huber lab scope, it is 283.385):
  mutate(nanoCONC = ((nanoAvg * 283.385)/(VOL * 0.01 * 0.9)),
         microCONC = ((microAvg * 283.385)/(VOL * 0.01 * 0.9)),
         eukCONC = ((euksAvg * 283.385)/(VOL * 0.01 * 0.9))
         ) %>%
  # left_join(expmeta) %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  data.frame
head(counts_cellsml_all)
##                SAMPLE    Site        Name    EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp  T0-Rep3        T0      Rep3  50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3       T15      Rep3  50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3       T20      Rep3  50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3       T40      Rep3  40
## 5       Piccard-Plume Piccard       Plume  T0-Rep1        T0      Rep1 150
## 6       Piccard-Plume Piccard       Plume  T0-Rep2        T0      Rep2 150
##   totalFOV   nanoAvg   nanoVar    nanoSd   microAvg   microVar   microSd
## 1       30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2       30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3       30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4       30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5       30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6       30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
##     euksAvg   euksVar   euksSd  nanoCONC microCONC   eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630   0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333  41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481   0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676   0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957  62.97444   6.99716  69.97160
## 6 0.2666667 0.2712644 1.041661  41.98296  13.99432  55.97728

Replicates belong to the same experiment for either Bag or IGT incubation. Below, modify these names and label new column with bag or igt. And create an average across replicates.

Average cells/ml across replicates, pivot to long format

counts_cellsml_avg <- counts_cellsml_all %>%
  select(Site, Name, TimePoint, Replicate, nanoCONC, microCONC, eukCONC) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>%
  select(-Replicate) %>%
  pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
  group_by(Site, Name, TimePoint, EXP_TYPE, IGT_REP, VARIABLE) %>%
  # Calculate mean, variance, SD, min, and max
  summarise(MEAN = mean(CONCENTRATION),
            VAR = var(CONCENTRATION),
            SD = sd(CONCENTRATION),
            SEM =(sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
            MIN = min(CONCENTRATION),
            MAX = max(CONCENTRATION),
            .groups = 'drop_last') %>%
  data.frame
head(counts_cellsml_avg)
##      Site        Name TimePoint EXP_TYPE IGT_REP  VARIABLE      MEAN VAR SD SEM
## 1 Piccard LotsOShrimp        T0      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 2 Piccard LotsOShrimp        T0      Bag     Bag microCONC   0.00000  NA NA  NA
## 3 Piccard LotsOShrimp        T0      Bag     Bag  nanoCONC 230.90630  NA NA  NA
## 4 Piccard LotsOShrimp       T15      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 5 Piccard LotsOShrimp       T15      Bag     Bag microCONC  41.98296  NA NA  NA
## 6 Piccard LotsOShrimp       T15      Bag     Bag  nanoCONC 188.92333  NA NA  NA
##         MIN       MAX
## 1 230.90630 230.90630
## 2   0.00000   0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5  41.98296  41.98296
## 6 188.92333 188.92333
# View(counts_cellsml_avg)
# NOTES on calculations:
# VAR = takes the sum of the squares of each value's deviation from the mean and divides by the number of such values minus one. This differs from the calculation of variance across an entire population in that the latter divides by the size of the dataset without subtracting one.
# SD = standard deviation of all values
# SEM = standard deviation of sampling distribution; standard deviation divided by the square root of the sample size.
# Euk cell counts, all and averaged
head(counts_cellsml_all)
##                SAMPLE    Site        Name    EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp  T0-Rep3        T0      Rep3  50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3       T15      Rep3  50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3       T20      Rep3  50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3       T40      Rep3  40
## 5       Piccard-Plume Piccard       Plume  T0-Rep1        T0      Rep1 150
## 6       Piccard-Plume Piccard       Plume  T0-Rep2        T0      Rep2 150
##   totalFOV   nanoAvg   nanoVar    nanoSd   microAvg   microVar   microSd
## 1       30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2       30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3       30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4       30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5       30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6       30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
##     euksAvg   euksVar   euksSd  nanoCONC microCONC   eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630   0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333  41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481   0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676   0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957  62.97444   6.99716  69.97160
## 6 0.2666667 0.2712644 1.041661  41.98296  13.99432  55.97728
head(counts_cellsml_avg)
##      Site        Name TimePoint EXP_TYPE IGT_REP  VARIABLE      MEAN VAR SD SEM
## 1 Piccard LotsOShrimp        T0      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 2 Piccard LotsOShrimp        T0      Bag     Bag microCONC   0.00000  NA NA  NA
## 3 Piccard LotsOShrimp        T0      Bag     Bag  nanoCONC 230.90630  NA NA  NA
## 4 Piccard LotsOShrimp       T15      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 5 Piccard LotsOShrimp       T15      Bag     Bag microCONC  41.98296  NA NA  NA
## 6 Piccard LotsOShrimp       T15      Bag     Bag  nanoCONC 188.92333  NA NA  NA
##         MIN       MAX
## 1 230.90630 230.90630
## 2   0.00000   0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5  41.98296  41.98296
## 6 188.92333 188.92333

Parse experiment type information

# Convert to long format and add column that reports IGT vs bag experiment
plot_euk_conc <- counts_cellsml_all %>%
  select(Site, Name, TimePoint, Replicate, ends_with("CONC")) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  pivot_longer(cols = ends_with("CONC"), names_to = "VARIABLE", values_to = "CONCENTRATION") %>%
  data.frame
# head(plot_euk_conc)
unique(plot_euk_conc$Name)
## [1] "LotsOShrimp"    "Plume"          "Shrimpocalypse" "BSW"           
## [5] "MustardStand"   "OMT"            "Rav2"           "ShrimpHole"    
## [9] "X18"
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
plot_euk_conc$SiteOrder <- factor(plot_euk_conc$Site, levels = site_ids, labels = site_fullname)
plot_euk_conc$NameOrder <- factor(plot_euk_conc$Name, levels = vent_ids, labels = vent_fullname)

Box plot to report all eukaryote cell counts.

# Code for plot
conc_boxplot <- ggplot(plot_euk_conc, aes(x = NameOrder, 
                                          y = CONCENTRATION, 
                                          group = NameOrder,
                                          fill = VARIABLE,
                                          shape = EXP_TYPE)) +
    geom_boxplot() + 
    # Do not color by time point
    geom_jitter(color = "black", size = 2, aes(fill = VARIABLE,
                                          shape = EXP_TYPE)) +
    scale_shape_manual(values = c(21,24)) +
    scale_fill_manual(values = c("#e7298a", "#fcbba1", "#c6dbef")) +
    coord_flip() +
    scale_y_log10() +
    # scale_y_log10(limits = c(10,1000), expand = c(0, 0)) +
    facet_grid(SiteOrder ~ EXP_TYPE, space = "free", scale = "free") +
    theme_bw() + 
  theme(axis.text.x = element_text(angle = 0, h = 1, vjust = 1),
        strip.background = element_blank(),
        legend.position = "right",
        legend.title = element_blank()) +
    labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
         title = "Distribution of all eukaryotic cell counts")

Eukaryote cell concentration (cells/ml) are lower in the background and plume samples compared to vent sites. ~300 cells/ml in background and plume compared to ~1000 cells per ml at the vent sites. These values are also consistent between each vent site (Von Damm and Piccard) and between Bag and IGT samples.

Boxplot represents the median (line in box) and the 1st and 3rd quartiles in the lower and upper hinges, respectively (25th and 75th percentiles). Black data points are outliers from the boxplot. Upper and lower whiskers represent the 1.5 * interquartile ranges. Pink data points are the values contributing to the boxplot (individial counts across replicates and time points.)

conc_boxplot
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis
## Warning: Removed 39 rows containing non-finite values (stat_boxplot).

eukCONC is the sum of micro and nano. Because there was a discrepency between the micro and nano cell counts, we plan to combine for most of the analysis. Here we show that the cell concentration across replicate samples was similar throughout experiments. And that the bag versus IGT experiment results were within range of one another.

Subset to plot from T0 only.

vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")

plot_euk_format <- plot_euk_conc %>%
  filter(TimePoint == "T0" & (VARIABLE == "eukCONC")) %>%
  group_by(SiteOrder, NameOrder, TimePoint, EXP_TYPE, VARIABLE) %>%
  summarise(avg_conc = mean(CONCENTRATION),
            SEM_conc = (sd(CONCENTRATION)/sqrt(length(CONCENTRATION))),
            .groups = "rowwise") %>%
  unite(EXPERIMENT, SiteOrder, NameOrder, EXP_TYPE, remove = FALSE) %>%
  data.frame

# Factor
plot_euk_format$Site_Order <- factor(plot_euk_format$SiteOrder, levels = site_fullname, labels = site_fullname)

# View(plot_euk_format)
euk_plot <- ggplot(plot_euk_format, aes(x = NameOrder, y = avg_conc, fill = Site_Order)) +
  geom_errorbar(aes(ymax = (avg_conc + SEM_conc), ymin = (avg_conc - SEM_conc)), width = 0.2) +
  geom_point(aes(fill = Site_Order), color = "black", stat = "identity", size = 3, shape = 23) +
  facet_grid(.~ Site_Order, space = "free", scales = "free") +
  scale_fill_manual(values = c("#1c9099", "#de2d26")) +
  theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank()) +
  labs(x = "", y = bquote("Eukaryote cells "~mL^-1),
       title = "")
euk_plot

3 Prokaryote cell concentration

Import counts from DAPI counts of bacteria and archaea.

prok <- read.delim("input-data/prokINSITU-counts-compiled.txt")
insitu_proks <- prok %>% 
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>% 
  group_by(SAMPLE, Site, Name) %>% 
  summarise(MEAN = mean(CELLML),
            SD = sd(CELLML),
            SEM = (sd(CELLML)/sqrt(length(CELLML))),
            .groups = "rowwise") %>% 
  data.frame
# head(insitu_proks)

3.1 Microscopy cells per ml (prok)

Plot in situ prokaryote values. Average across dives where applicable.

insitu_proks$Name_order <- factor(insitu_proks$Name, levels = c("BSW", "Plume", "Quakeplume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole", "HotChimlet1", "ShrimpGulley", "SouthofHotChimlet", "SouthofLungSnack", "ArrowLoop", "Bartizan", "Rav1"), labels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
insitu_proks$Site_order <- factor(insitu_proks$Site, levels = site_ids, labels = site_fullname)
# unique(insitu_proks$Name)
prok_plot <- ggplot(insitu_proks, aes(x = Name_order, y = MEAN)) +
  geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
  geom_point(stat = "identity", shape = 23, aes(fill = Site), size = 3) +
  facet_grid(.~ Site_order, space = "free", scales = "free") +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(y = bquote("Prokaryote cells "~mL^-1), x = "", title = "") +
  scale_y_log10() +
  theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank())
prok_plot

3.2 Comparison with previous MCR

Can include years 2012 and 2013, seeking clarification on these data. ## Prokaryote biomass-literature

Work on this! get FS xx vlaues from Julie

4 Determine grazing effect

Calculate FLP per eukaryotic cell over time. Goal is to make these calculations and then determine best fit line. Slope of best fit line is the grazing rate. Need to take into account euk cells with FLPs and then the euk cells withOUT FLPs, these will be zeroes to take into account for FLPs/euk averages.

Retain previously generated dataframes:

# Euk cell counts, all and averaged
head(counts_cellsml_all)
##                SAMPLE    Site        Name    EXPID TimePoint Replicate VOL
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp  T0-Rep3        T0      Rep3  50
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp T15-Rep3       T15      Rep3  50
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp T20-Rep3       T20      Rep3  50
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp T40-Rep3       T40      Rep3  40
## 5       Piccard-Plume Piccard       Plume  T0-Rep1        T0      Rep1 150
## 6       Piccard-Plume Piccard       Plume  T0-Rep2        T0      Rep2 150
##   totalFOV   nanoAvg   nanoVar    nanoSd   microAvg   microVar   microSd
## 1       30 0.3666667 0.3781609 1.2298958 0.00000000 0.00000000 0.0000000
## 2       30 0.3000000 0.2172414 0.9321832 0.06666667 0.06436782 0.5074163
## 3       30 0.3333333 0.2988506 1.0933445 0.00000000 0.00000000 0.0000000
## 4       30 0.1666667 0.1436782 0.7580980 0.00000000 0.00000000 0.0000000
## 5       30 0.3000000 0.2862069 1.0699662 0.03333333 0.03333333 0.3651484
## 6       30 0.2000000 0.1655172 0.8136762 0.06666667 0.06436782 0.5074163
##     euksAvg   euksVar   euksSd  nanoCONC microCONC   eukCONC
## 1 0.3666667 0.3781609 1.229896 230.90630   0.00000 230.90630
## 2 0.3666667 0.2402299 0.980265 188.92333  41.98296 230.90630
## 3 0.3333333 0.2988506 1.093345 209.91481   0.00000 209.91481
## 4 0.1666667 0.1436782 0.758098 131.19676   0.00000 131.19676
## 5 0.3333333 0.3678161 1.212957  62.97444   6.99716  69.97160
## 6 0.2666667 0.2712644 1.041661  41.98296  13.99432  55.97728
head(counts_cellsml_avg)
##      Site        Name TimePoint EXP_TYPE IGT_REP  VARIABLE      MEAN VAR SD SEM
## 1 Piccard LotsOShrimp        T0      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 2 Piccard LotsOShrimp        T0      Bag     Bag microCONC   0.00000  NA NA  NA
## 3 Piccard LotsOShrimp        T0      Bag     Bag  nanoCONC 230.90630  NA NA  NA
## 4 Piccard LotsOShrimp       T15      Bag     Bag   eukCONC 230.90630  NA NA  NA
## 5 Piccard LotsOShrimp       T15      Bag     Bag microCONC  41.98296  NA NA  NA
## 6 Piccard LotsOShrimp       T15      Bag     Bag  nanoCONC 188.92333  NA NA  NA
##         MIN       MAX
## 1 230.90630 230.90630
## 2   0.00000   0.00000
## 3 230.90630 230.90630
## 4 230.90630 230.90630
## 5  41.98296  41.98296
## 6 188.92333 188.92333
head(counts_occur)
##     DATE   SAMPLE   EXPID VOL MAG FOV nanoNoFLP microNoFLP nanoFLP microFLP
## 1 4/2/20 VD-Plume T0-Rep1 150 100   1         0          0       0        0
## 2 4/2/20 VD-Plume T0-Rep1 150 100   2         0          0       0        0
## 3 4/2/20 VD-Plume T0-Rep1 150 100   3         3          0       0        0
## 4 4/2/20 VD-Plume T0-Rep1 150 100   4         1          0       4        0
## 5 4/2/20 VD-Plume T0-Rep1 150 100   5         0          1       0        0
## 6 4/2/20 VD-Plume T0-Rep1 150 100   6         0          0       0        0
##             NOTES DateCompiled nanoFLP_occur microFLP_occur nanoTOTAL
## 1               0       6/7/20             0              0         0
## 2               0       6/7/20             0              0         0
## 3               0       6/7/20             0              0         3
## 4               0       6/7/20             1              0         2
## 5 elongated cell?       6/7/20             0              0         0
## 6               0       6/7/20             0              0         0
##   microTOTAL euksTOTAL
## 1          0         0
## 2          0         0
## 3          0         3
## 4          0         2
## 5          1         1
## 6          0         0

Supplementary plot showing the concentration of eukaryotic cells over time for each experiment.

# head(counts_cellsml_avg)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")

counts_cellsml_avg$SiteOrder <- factor(counts_cellsml_avg$Site, levels = site_ids, labels = site_fullname)
counts_cellsml_avg$NameOrder <- factor(counts_cellsml_avg$Name, levels = vent_ids, labels = vent_fullname)

# Plot trend line of euk cell count for all experiments
counts_cellsml_avg %>%
  filter(VARIABLE == "eukCONC") %>%
  unite("Experiment", NameOrder, IGT_REP, sep = "-", remove = FALSE) %>%
  ggplot(aes(x = TimePoint, y = MEAN, shape = EXP_TYPE, fill = NameOrder)) +
    geom_path(aes(group = Experiment)) +
    # geom_errorbar(aes(ymax = (MEAN + SD), ymin = (MEAN - SD)), width = 0.2) +
    geom_errorbar(aes(ymax = (MEAN + SEM), ymin = (MEAN - SEM)), width = 0.2) +
    geom_point(stat = "identity", size = 2, aes(shape = EXP_TYPE)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_brewer(palette = "Paired") +
    scale_y_log10() +
    facet_wrap(SiteOrder ~ EXP_TYPE, scales = "free") +
    theme_classic() + theme(strip.background = element_blank(), 
                            legend.title = element_blank(),
                            title = element_text(size = 7, face = "bold"),
                            axis.title = element_text(size = 9)) +
    labs(title = "Total euk cell counts for each experiment", y = bquote("Average eukaryote cells "~mL^-1), x = "Time point") +
  guides(fill=guide_legend(override.aes=list(shape=21)))

## Determine FLP per euk Isolate the euk cell counts with FLPs, these are comma separated and must be separated into rows. Using counts_occur data frame from above.

# Select nano and micro counts with FLPs
counts_sepflp <- counts_occur %>% 
  filter(!NOTES == "Discard") %>% 
  filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
  select(DATE, SAMPLE, EXPID, VOL, MAG, FOV, nanoFLP, microFLP) %>%
  separate_rows(microFLP, sep = ",", convert = TRUE) %>%
  separate_rows(nanoFLP, sep = ",", convert = TRUE) %>%
  replace_na(list(microFLP = 0, nanoFLP = 0)) %>% #Replace NAs with zero
  data.frame

## Check, see FOV 23, separated into rows.
# View(counts_sepflp %>%
  # filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))
# View(counts_occur %>%
       # filter(SAMPLE == "VD-Rav2" & EXPID == "T10-Rep1"))

Subset counts that are greater than 0, so only euk cells observed to have FLPs. Create new column that calculates FLP per Euk - by dividing by 1.

counts_flp <- counts_sepflp %>%
  select(SAMPLE, EXPID, nano_size = nanoFLP, micro_size = microFLP) %>%
  pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_FLP") %>%
  filter(num_of_FLP > 0) %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>%
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  summarise(total_FLP = sum(num_of_FLP),
            total_euks_wflp = n(),
            .groups = "rowwise") %>%
  data.frame

OUTPUT COLUMNS: (1) total_FLP = sum of FLPs found inside a euk cell (2) total_euks_wflp = number of euks counted with ingested FLP

Repeat above operation for euk cells without any FLP. Here, subset total number of observations where there was a euk cell without FLP. These need to be counted as euk cell without an FLP. > Below code repeats process and compiles with other FLP/euk cell data.

counts_flp_compiled <- counts_occur %>% 
  filter(!(NOTES == "Discard")) %>% #Discard bad counts
  filter(!(NOTES == "DTAF stain prevented counts of FLP, Euks only")) %>%
  type.convert(as.is = TRUE) %>% #modify str() for columns
  select(SAMPLE, EXPID, nano_size = nanoNoFLP, micro_size = microNoFLP) %>% #select non flp
  pivot_longer(cols = ends_with("_size"), names_to = "SizeFrac", values_to = "num_of_euks") %>%
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate"), sep = "-", remove = FALSE) %>%
  mutate(EXP_TYPE = case_when(
    grepl("IGT", Replicate) ~ "IGT",
    grepl("Rep", Replicate) ~ "Bag"
  )) %>%
  mutate(IGT_REP = case_when(
    EXP_TYPE == "IGT" ~ Replicate,
    EXP_TYPE == "Bag" ~ "Bag")) %>% 
  # filter(num_of_euks > 0) %>% # Remove observed zero counts
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  summarise(total_euks_noFLP = sum(num_of_euks),
            .groups = "rowwise") %>%
  # Join with FLP count information
  ## SAMPLE, EXPID, EXPTYPE, IGTREP, and SizeFrac variables should match
  left_join(counts_flp) %>% 
  replace_na(list(total_FLP = 0, total_euks_wflp = 0)) %>% #Replace NAs with zero
  data.frame
## Joining, by = c("SAMPLE", "EXPID", "EXP_TYPE", "IGT_REP", "SizeFrac")
## Check data frames
# View(counts_flp_compiled)
# View(counts_occur)
# (filter(counts_flp_compiled, SAMPLE == "VD-X18", EXPID == "T40-Rep2"))
# (filter(counts_flp_compiled, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1"))
# head(counts_flp_compiled)

Extract total euk value by adding across nano and micro, then combine back with nano and micro counts.

counts_flp_compiled_all <- counts_flp_compiled %>% 
  # Exclude size fraction:
  group_by(SAMPLE, EXPID, EXP_TYPE, IGT_REP) %>%
  summarise(total_euks_noFLP = sum(total_euks_noFLP),
            total_FLP = sum(total_FLP), 
            total_euks_wflp = sum(total_euks_wflp),
            .groups = "rowwise") %>% 
  add_column(SizeFrac = "total_euks") %>% #Add SizeFrac column
  bind_rows(counts_flp_compiled) %>% # Combine back with flp compiled list
  data.frame
# head(counts_flp_compiled_all)
# filter(counts_flp_compiled_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_compiled_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")

All three size fractions are represented, total, micro, and nano.

# unique(counts_flp_compiled_all$SizeFrac)
# unique(counts_flp_compiled_all$SAMPLE)
# head(counts_flp_compiled_all[1:2,])

4.1 Calculate FLP per euk cell calculation

Import metadata which has the extact minutes for each time point used.

metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
head(metadata)
##            SAMPLE TimePoint Minutes SiteOrigin REP
## 1 VD-MustardStand        T0       0       Vent Bag
## 2 VD-MustardStand       T10      10       Vent Bag
## 3 VD-MustardStand       T15      15       Vent Bag
## 4 VD-MustardStand       T20      20       Vent Bag
## 5 VD-MustardStand       T40      40       Vent Bag
## 6   Piccard-Plume        T0       0      Plume Bag

Add metadata information to FLP data, reformat sample namse, and calculate FLP per euk as the total FLP divided by the total number of euk cells counted.

counts_flp_calcs_all <- counts_flp_compiled_all %>% 
  # Add in metadata
  # IGTXb are replicate counts, include them as replicates!
  separate(EXPID, c("TimePoint", "REP"), sep = "-", remove = FALSE) %>% mutate(
    REP = ifelse(grepl("IGT5b", REP), "IGT5", REP),
    REP = ifelse(grepl("IGT4b", REP), "IGT4", REP),
    REP = ifelse(grepl("Bag", EXP_TYPE), "Bag", REP)) %>% 
  left_join(metadata, by = c("SAMPLE" = "SAMPLE", "TimePoint" = "TimePoint", "REP" = "REP")) %>% 
  separate(SAMPLE, c("Site", "Name"), sep = "-", remove = FALSE) %>%
  separate(EXPID, c("TimePoint", "Replicate_ID"), sep = "-", remove = FALSE) %>%
  ## Treat repeated IGT counts completely separate
  # group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, SizeFrac) %>%
  ## Treat repeated IGT counts as replicates (e.g., IGT4b and IGT4 == IGT4)
  group_by(SAMPLE, Site, Name, EXPID, TimePoint, Replicate_ID, EXP_TYPE, REP, SizeFrac) %>%
  # FLPperEuk is the total FLP divided by the total number of euk cells counted
  mutate(FLPperEuk = total_FLP/(sum(total_euks_noFLP, total_euks_wflp))) %>%
  unite("Experiment", Name, REP, sep = "-", remove = FALSE) %>%
  data.frame

# filter(counts_flp_calcs_all, SAMPLE == "VD-X18", EXPID == "T40-Rep2")
# filter(counts_flp_calcs_all, SAMPLE == "Piccard-Plume", EXPID == "T10-Rep1")

View output

# head(counts_flp_calcs_all[1:2,])
unique(counts_flp_calcs_all$Experiment)
##  [1] "LotsOShrimp-Bag"     "Plume-Bag"           "Shrimpocalypse-IGT3"
##  [4] "Shrimpocalypse-Bag"  "BSW-Bag"             "MustardStand-Bag"   
##  [7] "OMT-IGT4"            "Rav2-IGT4"           "Rav2-IGT5"          
## [10] "Rav2-Bag"            "ShrimpHole-Bag"      "X18-Bag"
# View(counts_flp_calcs_all)

COLS: Timepoint, Minutes = time point label, actual incubated minutes COLS: Replicate_ID, REP, and IGT_REP = full replicate identified for IGTs and Bags, designation of biological replicates, and designation of technical replicates for IGT experiments

4.2 Calculate linear regression to obtain slope

Use lm() function in R to calculate linear regression for each experiment. Slope equates to grazing rate. Function inputs the FLP per euk cell data, performs regression and then adds a column for slope and r-squared values.

calculate_lm <- function(df){
  regression_1 <- df %>%
  type.convert(as.is = TRUE) %>%
  ## Keep technical replicates separate for IGTs
  # group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
  # nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
  ## Combine technical replicates for IGTs
  group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
  nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
  mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
         tidied = map(lm_fit, tidy)) %>% 
  unnest(tidied) %>% 
  # select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, term, estimate) %>%
  select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  data.frame
#   colnames(regression_1) <- c("SAMPLE", "Site", "Experiment", "Name", "IGT_REP", 
# "SizeFrac", "INTERCEPT", "SLOPE")
  colnames(regression_1) <- c("SAMPLE", "Site", 
                              "Experiment", "Name", "REP",
                              "SizeFrac", "INTERCEPT", "SLOPE")
  out_regression <- df %>%
  # group_by(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac) %>%
  # nest(-SAMPLE, -Site, -Experiment, -Name, -IGT_REP, -SizeFrac) %>%
  group_by(SAMPLE, Site, Experiment, Name, REP, SizeFrac) %>%
  nest(-SAMPLE, -Site, -Experiment, -Name, -REP, -SizeFrac) %>%
  mutate(lm_fit = map(data, ~lm(FLPperEuk ~ Minutes, data = .)),
         glanced = map(lm_fit, glance)) %>% 
  unnest(glanced) %>% 
  # select(SAMPLE, Site, Experiment, Name, IGT_REP, SizeFrac, r.squared) %>% 
  select(SAMPLE, Site, Experiment, Name, REP, SizeFrac, r.squared) %>% 
  right_join(regression_1) %>% 
  right_join(df) %>% 
  data.frame
  out_regression$SITE <- factor(out_regression$Site, levels = c("VD", "Piccard"))
  out_regression$TYPE <- factor(out_regression$EXP_TYPE, levels = c("Bag", "IGT"))
  return(out_regression)
}

Note that an error may occur when running the below function. This is due to the fact that some experiments did not have replicates.

calcs_wslope_regression <- calculate_lm(counts_flp_calcs_all)
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk)`?
## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable

## Warning in summary.lm(x): essentially perfect fit: summary may be unreliable
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
# head(calcs_wslope_regression[1:2,])
# View(calcs_wslope_regression)

Use below commented out to recalculate one linear regression. Above function uses the nest() capability of tidyverse. Below, one experiment is subset to checek the value.

# Extract only plume-bag experiment from VD
# tmp_plume <- filter(counts_flp_calcs_all, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks")
# tmp_plume # View
# Perform linear regression
# lm_out <- lm(FLPperEuk ~ Minutes, data = tmp_plume)
# Check output
# summary(lm_out)
# lm_out$coefficients #Intercept=intercept #Minutes = SLOPE
# Compare with nested function output
# filter(calcs_wslope_regression, Experiment == "Plume-Bag") %>% filter(Site == "VD") %>% filter(SizeFrac == "total_euks") %>% head
# View(calcs_wslope_regression)
# head(calcs_wslope_regression)
# unique(calcs_wslope_regression$Site)

4.2.1 Plot linear regression trend

calcs_wslope_regression %>% 
  filter(SizeFrac == "total_euks") %>% 
  # filter(TYPE != "IGT") %>% 
  unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>% 
  ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
  geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
  geom_point(stat = "identity", color = "black", 
             size = 2, aes(shape = TYPE, fill = Site)) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
  facet_wrap(. ~ EXPERIMENT) +
  # Report r.squared
  geom_text(aes(x = 42, y = max(FLPperEuk), label = paste(round(SLOPE, 4))), 
            vjust = 1, hjust = 0, size = 3) +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(color = "black", size = 7),
                     legend.title = element_blank(),
                     legend.position = "right")

Data points represent the FLP per euk cells (based on total eukaryote cells counts). Y-axis represents the duration of incubation (in minutes). The dashed purple line reprents the slope and intercept of the experiment.

4.3 Remove Tf-IGT experiments

IGT experiment results appear to have bottle effect, especially in the final time point. Additionally, due to the lack of biological replicates in the IGT experiments, technical replicates are treated as biological replicates in the regression below.

IGT_lm_woTf <- counts_flp_calcs_all %>% 
  # Select only IGT experiments with total eukaryotes, remove Tf (T3)
  filter(SizeFrac == "total_euks") %>% 
  filter(EXP_TYPE == "IGT" & !(TimePoint == "T3")) %>% 
  add_column(IGT_cor = "rm Tf") %>% 
  data.frame

# Recalculate lm(), but keep 
igt_regression_noTf <- calculate_lm(IGT_lm_woTf) # Recalculate
## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk, 
##     IGT_cor)`?

## Warning: All elements of `...` must be named.
## Did you want `data = c(EXPID, TimePoint, Replicate_ID, EXP_TYPE, IGT_REP, total_euks_noFLP, 
##     total_FLP, total_euks_wflp, Minutes, SiteOrigin, FLPperEuk, 
##     IGT_cor)`?
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")
## Joining, by = c("SAMPLE", "Site", "Experiment", "Name", "REP", "SizeFrac")

Plot IGT grazing experiment with newly calculated grazing effect.

igt_regression_noTf %>% 
  # filter(SizeFrac == "total_euks") %>% 
  # filter(TYPE != "IGT") %>% 
  unite(EXPERIMENT, SITE, Experiment, sep = " ", remove = FALSE) %>% 
  ggplot(aes(x = Minutes, y = FLPperEuk, fill = Site, shape = TYPE, group = Experiment)) +
  geom_abline(aes(slope = SLOPE, intercept = INTERCEPT), color = "black", linetype = "dashed", size = 1) +
  geom_point(stat = "identity", color = "black", 
             size = 2, aes(shape = TYPE, fill = Site)) +
  scale_shape_manual(values = c(24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  labs(x = "Minutes", y = bquote("FLP"~eukaryote^-1), title = "Grazing experiment regression") +
  facet_wrap(. ~ EXPERIMENT) +
  # Report r.squared
  geom_text(aes(x = 5, y = max(FLPperEuk), label = paste(round(SLOPE, 4))), 
            vjust = 1, hjust = 0, size = 3) +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(color = "black", size = 7),
                     legend.title = element_blank(),
                     legend.position = "right")

4.4 Compile grazing experiment results

calcs_wslope_regression_update <- calcs_wslope_regression %>% 
  filter(TYPE != "IGT") %>% 
  bind_rows(igt_regression_noTf %>% select(-IGT_cor)) %>% 
  data.frame

# Factor
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
# Factor for shipboard
calcs_wslope_regression_update$SiteOrder <- factor(calcs_wslope_regression_update$Site, levels = site_ids, labels = site_fullname)
calcs_wslope_regression_update$NameOrder <- factor(calcs_wslope_regression_update$Name, levels = vent_ids, labels = vent_fullname)

# dim(calcs_wslope_regression); dim(calcs_wslope_regression_update)
head(calcs_wslope_regression_update)
##                SAMPLE    Site      Experiment        Name REP   SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
## 6       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
##    r.squared INTERCEPT        SLOPE    EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829  T0-Rep3        T0         Rep3      Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3       T15         Rep3      Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3       T20         Rep3      Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3       T40         Rep3      Bag
## 5 0.02810998 0.7435088  0.005362856  T0-Rep1        T0         Rep1      Bag
## 6 0.02810998 0.7435088  0.005362856  T0-Rep2        T0         Rep2      Bag
##   IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1     Bag                9         3               2       0       Vent
## 2     Bag                8         4               3      15       Vent
## 3     Bag                9         2               1      20       Vent
## 4     Bag                5         0               0      40       Vent
## 5     Bag                6         6               4       0      Plume
## 6     Bag                5         5               3       0      Plume
##   FLPperEuk    SITE TYPE SiteOrder      NameOrder
## 1 0.2727273 Piccard  Bag   Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard  Bag   Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard  Bag   Piccard          Plume
## 6 0.6250000 Piccard  Bag   Piccard          Plume

4.5 Check FLP control experiments

All incubations had control experiments run alongside them. This was to ensure added FLP did not decrease or change in concentration over time.

bac_ctrl <- read.delim("input-data/bac-counts-compiled.txt")
# dim(bac_ctrl)
dtaf <- bac_ctrl %>% 
  separate(SampleID, c("exp", "Replicate", "TimePoint"), sep = "-", remove = FALSE) %>% 
  separate(Site, c("Site", "Name"), sep = "-", remove = FALSE) %>% 
  filter(Stain == "DTAF") %>% 
  data.frame
## Warning: Expected 2 pieces. Additional pieces discarded in 17 rows [33, 34, 35,
## 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49].
# View(bac_ctrl)
# head(dtaf)

dtaf_avg <- dtaf %>% 
  group_by(TimePoint, Stain, Site, Name) %>% 
  summarise(Avg_cellsperml = mean(Cells.ml)) %>% 
  data.frame
## `summarise()` has grouped output by 'TimePoint', 'Stain', 'Site'. You can override using the `.groups` argument.
dtaf_avg %>% 
  filter(Site != "IGT") %>% 
  ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
  geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site != "IGT"), aes(
                                           ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
                                           ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
  geom_line(aes(group = Name)) +
  geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
  # scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
  facet_wrap(Name ~ Site) +
  scale_y_log10() +
  theme_bw() + theme(strip.background = element_blank(), 
                          legend.title = element_blank(),
                     axis.text = element_text(size = 10, color = "black"),
                          title = element_text(size = 10, face = "bold"),
                          axis.title = element_text(size = 9)) +
  labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")

dtaf_avg %>% 
  filter(Site == "IGT") %>% 
  ggplot(aes(x = TimePoint, y = Avg_cellsperml, fill = Name, shape = Site)) +
  geom_rect(data = filter(dtaf_avg, TimePoint == "T0", Site == "IGT"), aes(
                                           ymin = (Avg_cellsperml-(0.1*Avg_cellsperml)),
                                           ymax = (Avg_cellsperml+(0.1*Avg_cellsperml))), color = NA, alpha = 0.4, xmin = 0, xmax = 6, fill = "black") +
  geom_line(aes(group = Name)) +
  geom_point(stat = "identity", aes(shape = Site, fill = Name), size = 2) +
  # scale_fill_manual(values = c("black","#9970ab", "#5aae61")) +
  facet_wrap(Name ~ Site) +
  scale_y_log10() +
  theme_bw() + theme(strip.background = element_blank(), 
                          legend.title = element_blank(),
                     axis.text = element_text(size = 10, color = "black"),
                          title = element_text(size = 10, face = "bold"),
                          axis.title = element_text(size = 9)) +
  labs(title = "FLP counts for controls", y = bquote("Log FLP "~mL^-1), x = "Time point")

4.6 Table with grazing values

head(calcs_wslope_regression_update)
##                SAMPLE    Site      Experiment        Name REP   SizeFrac
## 1 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 2 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 3 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 4 Piccard-LotsOShrimp Piccard LotsOShrimp-Bag LotsOShrimp Bag total_euks
## 5       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
## 6       Piccard-Plume Piccard       Plume-Bag       Plume Bag total_euks
##    r.squared INTERCEPT        SLOPE    EXPID TimePoint Replicate_ID EXP_TYPE
## 1 0.66025258 0.3517002 -0.007605829  T0-Rep3        T0         Rep3      Bag
## 2 0.66025258 0.3517002 -0.007605829 T15-Rep3       T15         Rep3      Bag
## 3 0.66025258 0.3517002 -0.007605829 T20-Rep3       T20         Rep3      Bag
## 4 0.66025258 0.3517002 -0.007605829 T40-Rep3       T40         Rep3      Bag
## 5 0.02810998 0.7435088  0.005362856  T0-Rep1        T0         Rep1      Bag
## 6 0.02810998 0.7435088  0.005362856  T0-Rep2        T0         Rep2      Bag
##   IGT_REP total_euks_noFLP total_FLP total_euks_wflp Minutes SiteOrigin
## 1     Bag                9         3               2       0       Vent
## 2     Bag                8         4               3      15       Vent
## 3     Bag                9         2               1      20       Vent
## 4     Bag                5         0               0      40       Vent
## 5     Bag                6         6               4       0      Plume
## 6     Bag                5         5               3       0      Plume
##   FLPperEuk    SITE TYPE SiteOrder      NameOrder
## 1 0.2727273 Piccard  Bag   Piccard Lots 'O Shrimp
## 2 0.3636364 Piccard  Bag   Piccard Lots 'O Shrimp
## 3 0.2000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 4 0.0000000 Piccard  Bag   Piccard Lots 'O Shrimp
## 5 0.6000000 Piccard  Bag   Piccard          Plume
## 6 0.6250000 Piccard  Bag   Piccard          Plume
# Generate final table
table_grazerate <- calcs_wslope_regression_update %>% 
  filter(SizeFrac == "total_euks") %>% 
  group_by(SAMPLE, SiteOrder, NameOrder, SLOPE, EXP_TYPE) %>% 
  mutate(No_of_Replicates = n_distinct(Replicate_ID)) %>%
  select(SAMPLE, Site=SiteOrder, Name=NameOrder, RATE = SLOPE, EXP_TYPE, No_of_Replicates) %>% 
  distinct() %>% 
  data.frame
# View(table_grazerate)
# write_delim(table_grazerate, path = "Table-grazingrates.txt", delim = "\t")

Add FLP conc to table

# head(table_grazerate)
dtaf_igt <- 5352.8278 # Manually insert FLP concentration for IGT experiments; this value is estimated from how IGT FLP spike-ins were calculated

table_grazerate_wflp <- bac_ctrl %>% 
  filter(FLP_t0 == "use") %>% 
  add_column(EXP_TYPE = "Bag") %>% 
  group_by(Site, EXP_TYPE) %>% 
  summarise(FLP_conc = mean(Cells.ml)) %>% 
  right_join(table_grazerate, by = c("Site" = "SAMPLE", "EXP_TYPE" = "EXP_TYPE")) %>% 
  mutate(FLP_conc = ifelse(EXP_TYPE == "IGT", dtaf_igt, FLP_conc)) %>% 
  data.frame
## `summarise()` has grouped output by 'Site'. You can override using the `.groups` argument.
# View(table_grazerate_wflp)

Negative grazing rates equal undetected (change to zero).

# head(table_grazerate)
bsw <- c("Plume", "Background")
grazerate_plot <- table_grazerate %>% 
  type.convert(as.is = TRUE) %>% 
  mutate(detected = case_when(
    RATE < 0 ~ "Not detected",
    TRUE ~ "Detected")) %>% 
  mutate(type = case_when(
    Name %in% bsw ~ Name,
    TRUE ~ paste("Vent", EXP_TYPE, sep="-")
  )) %>% 
  mutate(GRAZE_RATE = case_when(
    RATE < 0 ~ 0,
    TRUE ~ RATE
  )) %>% 
  mutate(type_site = case_when(
    Name %in% bsw ~ Name,
    TRUE ~ "Vent"
  )) %>% 
  # filter(detected == "Detected") %>% 
  data.frame
# Factoring
type_order <- c("Vent-Bag", "Vent-IGT", "Plume", "Background")
grazerate_plot$TYPE <- factor(grazerate_plot$type, levels = type_order)
vent_ids <- c("BSW","Plume", "Shrimpocalypse", "LotsOShrimp", "X18", "OMT", "Rav2", "MustardStand", "ShrimpHole")
vent_fullname <- c("Background","Plume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole")
site_ids <- c("VD", "Piccard")
site_fullname <- c("Von Damm", "Piccard")
vent_colors <- c("#0868ac", "#41ab5d", "#e7298a", "#c994c7", "#fc4e2a", "#fed976", "#6a51a3", "#ffeda0", "#a1d99b")
names(vent_colors) <- vent_fullname
grazerate_plot$NAME <- factor(grazerate_plot$Name, levels = vent_fullname)
# svg("", h =, w = )
grazing_min_plot <- grazerate_plot %>% 
  ggplot(aes(y = GRAZE_RATE, x = NAME, shape = EXP_TYPE, fill = Site)) +
  geom_jitter(stat = "identity", aes(shape = EXP_TYPE, fill = Site),
              color = "black", size = 3, width = 0.3) +
  scale_shape_manual(values = c(21, 24)) +
  scale_fill_manual(values = c("#de2d26", "#1c9099")) +
  facet_grid(.~Site, space = "free", scales = "free") +
  # coord_flip() +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
           panel.background = element_blank(), 
           axis.line = element_line(colour = "black"), 
           axis.text.x = element_text(color="black", size = 12, 
                                      angle = 45, hjust = 1, vjust = 1), 
           axis.text.y = element_text(color="black", size = 12),
           axis.title =element_text(color="black", size = 12),
           axis.ticks = element_line(),
           strip.text =element_blank(), legend.title = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
       shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = bquote("FLP " ~grazer^-1 ~min^-1))
# dev.off()
grazing_min_plot

5 Plot cell conc & grazing

Generate final plots that show eukaryote and prokaryote biomass and cell concentration. Save final table reporting carbon estimates.

# Subset the average in situ prok cells/ml for non-background samples
tmp <- filter(insitu_proks, Name != "BSW", Name != "Plume") %>% select(MEAN)
avg_insitu <- mean(tmp$MEAN)
# head(insitu_proks)
# Add to master table with data
table_grazerate_wflp_wprok <- insitu_proks %>% 
  select(Site = SAMPLE, Prok_conc = MEAN, Prok_sem = SEM) %>% 
  right_join(table_grazerate_wflp) %>% 
  mutate(Prok_conc = ifelse(is.na(Prok_conc), avg_insitu, Prok_conc)) %>% 
  data.frame
## Joining, by = "Site"

Table with all cell count data

# head(table_grazerate_wflp_wprok)
table_grazerate_wflp_wprok_weuk <- plot_euk_format %>% 
  select(Name = NameOrder, Site.y.y = SiteOrder, euk_conc = avg_conc, EXP_TYPE, euk_conc_sem = SEM_conc) %>% 
  right_join(table_grazerate_wflp_wprok) %>% 
  select(SITE = Site.y.y, NAME = Name, EXP = EXP_TYPE, SAMPLE = Site, RATE_min = RATE, REPs = No_of_Replicates, FLP_ml = FLP_conc, PROK_ml = Prok_conc, PROK_sem = Prok_sem, EUK_ml = euk_conc, EUK_sem = euk_conc_sem) %>% 
  mutate(RATE_min = case_when(
    RATE_min < 0 ~ 0,
    TRUE ~ RATE_min
  )) %>% 
  data.frame
## Joining, by = c("Name", "Site.y.y", "EXP_TYPE")

6 Grazing rate calculations

head(table_grazerate_wflp_wprok_weuk)
##       SITE         NAME EXP   SAMPLE    RATE_min REPs    FLP_ml   PROK_ml
## 1 Von Damm   Background Bag   VD-BSW 0.002958889    3  4861.824  37889.62
## 2 Von Damm        Plume Bag VD-Plume 0.005274231    3 34242.354  16478.31
## 3 Von Damm         X-18 Bag   VD-X18 0.001744429    2  3148.722 111429.78
## 4 Von Damm Old Man Tree IGT   VD-OMT 0.014510943    2  5352.828  71147.40
## 5 Von Damm   Ravelin #2 Bag  VD-Rav2 0.003470217    3 51890.942  71147.40
## 6 Von Damm   Ravelin #2 IGT  VD-Rav2 0.015395240    2  5352.828  71147.40
##   PROK_sem    EUK_ml   EUK_sem
## 1 8608.427  91.83773  21.86613
## 2 2623.935 157.77468  67.09859
## 3 2973.793 314.87222 104.95741
## 4       NA 472.30833 122.45031
## 5       NA 409.33389  73.47019
## 6       NA 620.99799 123.17702

Based on Unrein et al. 2007, we use the estimated grazing rate, in situ prok abundance, in situ euk abundance, and the concentration of FLP to make additional estimates.

table_wcalcs <- table_grazerate_wflp_wprok_weuk %>%
  # Ingestion rate per hour
  mutate(RATE_hr = (RATE_min * 60),
         RATE_day = (RATE_hr * 24), #Compare to GR?
         # FLP concentration per L
         FLP_L = (FLP_ml * 1000),
         # mL per grazer per hr
         CLEARANCE_RATE_ml = (RATE_hr/FLP_ml),
         # nL per grazer per hour
         CLEARANCE_RATE_nL = ((RATE_hr/FLP_ml)/1.00E+6), 
         # proks per grazer per hr
         SPEC_GRAZE_RATE_hr = (CLEARANCE_RATE_ml * PROK_ml), 
         # proks per grazer per day
         GRAZE_RATE_DAY = (24 * SPEC_GRAZE_RATE_hr),
         # proks per ml per hr
         GRAZING_EFFECT_hr = (SPEC_GRAZE_RATE_hr * EUK_ml),
         GRAZING_EFFECT_hr_min = (SPEC_GRAZE_RATE_hr * (EUK_ml - EUK_sem)),
         GRAZING_EFFECT_hr_max = (SPEC_GRAZE_RATE_hr * (EUK_ml + EUK_sem)),
         # cells per ml per day
         GRAZING_EFFECT_day = ((SPEC_GRAZE_RATE_hr * 24) * EUK_ml),
         # Percentage per day
         BAC_TURNOVER_PERC = 100*(GRAZING_EFFECT_day / PROK_ml),
         BAC_TURNOVER_PERC_min = 100*(GRAZING_EFFECT_day / (PROK_ml - PROK_sem)),
         BAC_TURNOVER_PERC_max = 100*(GRAZING_EFFECT_day / (PROK_ml + PROK_sem))) %>% 
  data.frame
# View(table_wcalcs)

Explaination of units for table with calculated values. * RATE_min & RATE_hr = Grazing rate as ‘FLPs per grazer per minute’ and per hour * CLEARANCE_RATE = ml or nL per grazer per hour * SPEC_GRAZE_RATE (Specific grazing rate) = Prokaryotes per grazer per hour * GRAZING EFFECT = bacteria per ml per hour * Bacterial turnover rate = % per day # Carbon biomass table construction

Unrein, F., Massana, R., Alonso-Sáez, L., and Gasol, J.M. (2007) Significant year-round effect of small mixotrophic flagellates on bacterioplankton in an oligotrophic coastal system. Limnol Oceanogr 52: 456–469.

7 Estimates of carbon biomass and consumption

7.1 Protist carbon-measured

References for estimating biovolume Pernice, M.C., Forn, I., Gomes, A., Lara, E., Alonso-Sáez, L., Arrieta, J.M., et al. (2015) Global abundance of planktonic heterotrophic protists in the deep ocean. ISME J 9: 782–792.

# Import manual biovolume measurements
biov <- read.delim("input-data/biovol-euk-12-10-2020.txt")
# head(biov)

# Calculate volume
biov_calc <- biov %>% 
  mutate(SizeFrac = case_when(
    h >= 20 ~ "micro",
    TRUE ~ "nano")) %>% 
  mutate(Volume = ((pi/6) * (d^2) * d)) %>% # Calculate volume (um cubed) # Hillebrand et al. 1999
  mutate(pgC_cell = (0.216 * (Volume^0.939))) %>% # Calculate Cell biomass in pg C per cell # Menden-Deuer and Lessard 2000
  data.frame
# View(biov_calc)
biov_calc
##    EXP VENT_BSW      h      d SizeFrac     Volume     pgC_cell
## 1  IGT     vent 30.077 25.764    micro 8954.44130 1110.2426245
## 2  IGT     vent 89.582 10.000    micro  523.59878   77.1957618
## 3  Bag      BSW 14.595  8.036     nano  271.71800   41.6956679
## 4  Bag      BSW 12.480  8.982     nano  379.41786   57.0486460
## 5  Bag     vent  9.218  3.120     nano   15.90239    2.9015292
## 6  IGT     vent 17.255  9.986     nano  521.40274   76.8917043
## 7  Bag     vent 41.153 21.000    micro 4849.04826  624.1445904
## 8  IGT     vent 10.282  4.136     nano   37.04591    6.4194942
## 9  IGT     vent 29.776 25.852    micro 9046.50993 1120.9583343
## 10 IGT     vent 10.991  4.000     nano   33.51032    5.8424695
## 11 Bag     vent 14.333  2.000     nano    4.18879    0.8290772
## 12 Bag     vent 36.164  3.000    micro   14.13717    2.5980292
## 13 Bag      BSW 16.206 14.924     nano 1740.42111  238.4669404
## 14 Bag      BSW  7.000  7.000     nano  179.59438   28.2640658
## 15 Bag     vent 10.069  5.000     nano   65.44985   10.9544849

Volume is reported as um^3

# Volume by experiment type
biov_calc %>% 
  group_by(EXP) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 x 3
##   EXP     VOL     C
## * <fct> <dbl> <dbl>
## 1 Bag    836.  112.
## 2 IGT   3186.  400.
# Volume by euk size
biov_calc %>% 
  group_by(SizeFrac) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 x 3
##   SizeFrac   VOL     C
## * <chr>    <dbl> <dbl>
## 1 micro    4678. 587. 
## 2 nano      325.  46.9
# Volume by site
biov_calc %>% 
  group_by(VENT_BSW) %>% summarise(VOL = mean(Volume), C = mean(pgC_cell))
## # A tibble: 2 x 3
##   VENT_BSW   VOL     C
## * <fct>    <dbl> <dbl>
## 1 BSW       643.  91.4
## 2 vent     2188. 276.
# head(biov_calc)
euk_vol <- mean(biov_calc$Volume);euk_vol # in um^3
## [1] 1775.759
euk_carbon <- mean(biov_calc$pgC_cell); euk_carbon # in pg C per cell
## [1] 226.9636
euk_carbon_min <- min(biov_calc$pgC_cell); euk_carbon_min
## [1] 0.8290772
euk_carbon_max <- max(biov_calc$pgC_cell); euk_carbon_max
## [1] 1120.958
# euk_carbon

Avg euk biomass pg C per individual cell == {r}euk_carbon

7.2 Protist carbon-literature

Compare with Menden-Deuer and Lessard 2000, Table 2 - using only the heterotrophic species measured. Based on Table 2, the min volume was 4745 and the maximum was 1.2 x10^7 µm^3. Carbon content was measured at pg per cell, this was 469.48-35,339 pg per cell.

Import the heterotroph species volume and carbon content to compare to my measured values.

# Hu-measured
range(biov_calc$Volume)
## [1]    4.18879 9046.50993
range(biov_calc$pgC_cell)
## [1]    0.8290772 1120.9583343
c_prev <- read.delim("input-data/md-lessard-2000.txt") # Table 2, heterotrophs only

# combine and plot
carbon_compare <- c_prev %>% 
  add_column(source = "Menden-Deuer Lessard") %>% 
  select(source, Volume = vol, pgC_cell) %>% 
  rbind(biov_calc %>% add_column(source = "MCR") %>% select(source, Volume, pgC_cell)) %>% 
    ggplot(aes(x = Volume, y = pgC_cell, fill = source)) +
      geom_point(aes(fill = source), shape = 23, color = "black", size = 3) +
      scale_y_log10() + scale_x_log10() +
      labs(title = "Compare literature to measured cell volume & C content",
           x = bquote("Volume" ~µm^-3),
           y = bquote("pg C" ~cell^-1)) +
      theme_bw() + theme(legend.title = element_blank(),
                         axis.title = element_text(size = 14),
                         axis.text = element_text(size = 14),
                         legend.text = element_text(size = 14))

carbon_compare

Upon comparison, the measured carbon content was much lower from the grazing experiments. This makes sense, as I am looking at preserved specimen and a smaller total number of cells. AND the deep-sea protist cell sizes may be smaller overall.

Find lowest estimates or protist carbon, benthic estimates, and others? How does it compare to my measurements?

euk_carbon_lit_mean <- mean(c_prev$pgC_cell)
euk_carbon_lit_min <- min(c_prev$pgC_cell)
euk_carbon_lit_max <- max(c_prev$pgC_cell)

7.3 Prokaryote carbon-literature

Below adding in biomass estimates from prokaryotes and protists.

bac_carbon_ug <- (86)*(1.00E-9) # From Derived from Morono et al. 2011 
# bac_carbon_ug
bac_carbon_ug_2 <- (173)*(1.00E-9) # Derived from McNichol et al. 2018; LOFERER-KRO ̈ ßBACHER, J. KLIMA & R. PSENNER 1998
# table_wcalcs

7.4 Generate carbon biomass table

Incorporate calculations that include biomass of population and ug C consumed. For rate measurements, only incorporate the Morono et al. 2011 biomass for prokaryotes. This way it is on the lower end and is comparable to Gorda Ridge work.

NEED TO MODIFY SO IT INCLUDES BIOMASS OF EUK AND PROK AT EACH SITE

bsw <- c("Plume", "Background")

table_wcalcs_biomass <- table_wcalcs %>% 
  add_column(euk_C_ug_Hu = (euk_carbon / (1.00E+06))) %>% # Convert to ug from pg
  add_column(euk_C_ug_lit = (euk_carbon_lit_mean / (1.00E+06))) %>% # literature
  add_column(bac_C_ug = bac_carbon_ug) %>% 
  add_column(bac_C_ug_2 = bac_carbon_ug_2) %>%
  # Grazing rate in ug C per bac per day
  mutate(RATE_ugCbac_pergrazer_perday = (RATE_hr * 24 * bac_C_ug), # Grazing rate as ug C per grazer per day
         # % of cell carbon per day
         SPEC_INGESTION_RATE = (RATE_ugCbac_pergrazer_perday / euk_C_ug_Hu),
         SPEC_INGESTION_RATE_lit = (RATE_ugCbac_pergrazer_perday / euk_C_ug_lit),
         Prok_biomass = PROK_ml * bac_carbon_ug,
         Euk_biomass_Hu = EUK_ml * euk_C_ug_Hu,
         Euk_biomass_lit = EUK_ml * euk_C_ug_lit,
         Prok_biomass_L = PROK_ml * bac_carbon_ug * 1000,
         Euk_biomass_Hu_L = EUK_ml * euk_C_ug_Hu * 1000,
         Euk_biomass_lit_L = EUK_ml * euk_C_ug_lit * 1000,
         # Repeat with SEM values
         Prok_biomass_sem = PROK_sem * bac_carbon_ug,
         Euk_biomass_Hu_sem = EUK_sem * euk_C_ug_Hu,
         Euk_biomass_lit_sem = EUK_sem * euk_C_ug_lit,
         Prok_biomass_sem_L = PROK_sem * (bac_carbon_ug* 1000),
         Euk_biomass_Hu_sem_L = EUK_sem * (euk_C_ug_Hu * 1000),
         Euk_biomass_lit_sem_L = EUK_sem * (euk_C_ug_lit * 1000)) %>% 
  type.convert(as.is = TRUE) %>%
  mutate(detected = case_when(
    RATE_min < 0 ~ "Not detected",
    TRUE ~ "Detected")) %>% 
  mutate(type = case_when(
    NAME %in% bsw ~ NAME,
    TRUE ~ paste("Vent", EXP, sep="-")
  )) %>% 
  mutate(GRAZE_RATE = case_when(
    RATE_min < 0 ~ 0,
    TRUE ~ RATE_min
  )) %>% 
  mutate(type_site = case_when(
    NAME %in% bsw ~ NAME,
    TRUE ~ "Vent"
  )) %>%
  data.frame
# View(table_wcalcs_biomass)
  • Grazing rate column == FLP per minute consumed
  • Grazing effect hr == cells per ml per hr

Also make a “bounded” table that demonstrates the ug C consumed in the context of McNichol et al.

# G = number of cells grazed during experiment duration
table_wcalcs_biomass_bounded <- table_wcalcs_biomass %>% 
  add_column(fgC_cell = 86) %>% # Add in Morono et al. 2011 value
  mutate(
    # cells_consumed_perday = (G / 1), # Rate of cells consumed * in situ prok, per day
    fgC_ml_perday = (GRAZING_EFFECT_day * fgC_cell), # Convert cell amount to fg C
    ugC_L_perday = (fgC_ml_perday * (1e-09) * 1000), # Convert to ug C per L
    lower_mcnichol = 100*(ugC_L_perday / 17.3),
    upper_mcnichol = 100*(ugC_L_perday / 321.4)
  ) %>% 
  data.frame
head(table_wcalcs_biomass_bounded)
##       SITE         NAME EXP   SAMPLE    RATE_min REPs    FLP_ml   PROK_ml
## 1 Von Damm   Background Bag   VD-BSW 0.002958889    3  4861.824  37889.62
## 2 Von Damm        Plume Bag VD-Plume 0.005274231    3 34242.354  16478.31
## 3 Von Damm         X-18 Bag   VD-X18 0.001744429    2  3148.722 111429.78
## 4 Von Damm Old Man Tree IGT   VD-OMT 0.014510943    2  5352.828  71147.40
## 5 Von Damm   Ravelin #2 Bag  VD-Rav2 0.003470217    3 51890.942  71147.40
## 6 Von Damm   Ravelin #2 IGT  VD-Rav2 0.015395240    2  5352.828  71147.40
##   PROK_sem    EUK_ml   EUK_sem   RATE_hr  RATE_day    FLP_L CLEARANCE_RATE_ml
## 1 8608.427  91.83773  21.86613 0.1775333  4.260800  4861824      3.651579e-05
## 2 2623.935 157.77468  67.09859 0.3164538  7.594892 34242354      9.241591e-06
## 3 2973.793 314.87222 104.95741 0.1046657  2.511977  3148722      3.324070e-05
## 4       NA 472.30833 122.45031 0.8706566 20.895758  5352828      1.626536e-04
## 5       NA 409.33389  73.47019 0.2082130  4.997113 51890942      4.012512e-06
## 6       NA 620.99799 123.17702 0.9237144 22.169145  5352828      1.725657e-04
##   CLEARANCE_RATE_nL SPEC_GRAZE_RATE_hr GRAZE_RATE_DAY GRAZING_EFFECT_hr
## 1      3.651579e-11          1.3835695      33.205668         127.06388
## 2      9.241591e-12          0.1522858       3.654860          24.02685
## 3      3.324070e-11          3.7040034      88.896081        1166.28778
## 4      1.626536e-10         11.5723785     277.737084        5465.73080
## 5      4.012512e-12          0.2854798       6.851515         116.85656
## 6      1.725657e-10         12.2775993     294.662382        7624.36451
##   GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max GRAZING_EFFECT_day
## 1              96.81058             157.31719          3049.5332
## 2              13.80868              34.24501           576.6444
## 3             777.52519            1555.05037         27990.9067
## 4            4048.68948            6882.77212        131177.5392
## 5              95.88230             137.83081          2804.5574
## 6            6112.04639            9136.68264        182984.7483
##   BAC_TURNOVER_PERC BAC_TURNOVER_PERC_min BAC_TURNOVER_PERC_max  euk_C_ug_Hu
## 1          8.048465             10.414647              6.558411 0.0002269636
## 2          3.499414              4.162182              3.018725 0.0002269636
## 3         25.119772             25.808540             24.466811 0.0002269636
## 4        184.374334                    NA                    NA 0.0002269636
## 5          3.941897                    NA                    NA 0.0002269636
## 6        257.191065                    NA                    NA 0.0002269636
##   euk_C_ug_lit bac_C_ug bac_C_ug_2 RATE_ugCbac_pergrazer_perday
## 1   0.01147298  8.6e-08   1.73e-07                 3.664288e-07
## 2   0.01147298  8.6e-08   1.73e-07                 6.531607e-07
## 3   0.01147298  8.6e-08   1.73e-07                 2.160300e-07
## 4   0.01147298  8.6e-08   1.73e-07                 1.797035e-06
## 5   0.01147298  8.6e-08   1.73e-07                 4.297517e-07
## 6   0.01147298  8.6e-08   1.73e-07                 1.906547e-06
##   SPEC_INGESTION_RATE SPEC_INGESTION_RATE_lit Prok_biomass Euk_biomass_Hu
## 1         0.001614483            3.193840e-05  0.003258508     0.02084382
## 2         0.002877822            5.693033e-05  0.001417135     0.03580910
## 3         0.000951827            1.882946e-05  0.009582961     0.07146452
## 4         0.007917726            1.566319e-04  0.006118676     0.10719678
## 5         0.001893483            3.745771e-05  0.006118676     0.09290388
## 6         0.008400232            1.661770e-04  0.006118676     0.14094392
##   Euk_biomass_lit Prok_biomass_L Euk_biomass_Hu_L Euk_biomass_lit_L
## 1        1.053653       3.258508         20.84382          1053.653
## 2        1.810146       1.417135         35.80910          1810.146
## 3        3.612524       9.582961         71.46452          3612.524
## 4        5.418786       6.118676        107.19678          5418.786
## 5        4.696281       6.118676         92.90388          4696.281
## 6        7.124700       6.118676        140.94392          7124.700
##   Prok_biomass_sem Euk_biomass_Hu_sem Euk_biomass_lit_sem Prok_biomass_sem_L
## 1     0.0007403247        0.004962814           0.2508697          0.7403247
## 2     0.0002256584        0.015228935           0.7698210          0.2256584
## 3     0.0002557462        0.023821507           1.2041746          0.2557462
## 4               NA        0.027791758           1.4048704                 NA
## 5               NA        0.016675055           0.8429222                 NA
## 6               NA        0.027956696           1.4132080                 NA
##   Euk_biomass_Hu_sem_L Euk_biomass_lit_sem_L detected       type  GRAZE_RATE
## 1             4.962814              250.8697 Detected Background 0.002958889
## 2            15.228935              769.8210 Detected      Plume 0.005274231
## 3            23.821507             1204.1746 Detected   Vent-Bag 0.001744429
## 4            27.791758             1404.8704 Detected   Vent-IGT 0.014510943
## 5            16.675055              842.9222 Detected   Vent-Bag 0.003470217
## 6            27.956696             1413.2080 Detected   Vent-IGT 0.015395240
##    type_site fgC_cell fgC_ml_perday ugC_L_perday lower_mcnichol upper_mcnichol
## 1 Background       86     262259.86   0.26225986      1.5159529     0.08159921
## 2      Plume       86      49591.42   0.04959142      0.2866556     0.01542981
## 3       Vent       86    2407217.98   2.40721798     13.9145548     0.74897884
## 4       Vent       86   11281268.37  11.28126837     65.2096438     3.51003994
## 5       Vent       86     241191.93   0.24119193      1.3941730     0.07504416
## 6       Vent       86   15736688.36  15.73668836     90.9635165     4.89629383
# View(table_wcalcs_biomass_bounded)
write_delim(table_wcalcs_biomass_bounded, path = "table-wcalc.txt", delim = "\t")
# ?write_delim()

7.5 Compilation of plots

Format and factor values to plot. * Grazing rate column == FLP per minute consumed * Grazing effect hr == cells per ml per hr * Specific ingestion rate == % of cell carbon per day, estimated with literature and measured carbon values

biomass_rate_plot <- table_wcalcs_biomass_bounded %>% 
  select(SITE, NAME, EXP, SAMPLE, type, 
        PROK_ml, EUK_ml, PROK_sem, EUK_sem,
        Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
        Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
        GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, BAC_TURNOVER_PERC_min, BAC_TURNOVER_PERC_max,
        GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
        SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit) %>% 
  pivot_longer(cols = c(PROK_ml, EUK_ml,
                        Prok_biomass_L, Euk_biomass_Hu_L, Euk_biomass_lit_L,
                        GRAZING_EFFECT_hr, BAC_TURNOVER_PERC,
                        GRAZE_RATE, RATE_ugCbac_pergrazer_perday,
                        SPEC_INGESTION_RATE, SPEC_INGESTION_RATE_lit),
               names_to = "Variable", values_to = "Value") %>% 
  pivot_longer(cols = c(PROK_sem, EUK_sem, Prok_biomass_sem_L, Euk_biomass_Hu_sem_L, Euk_biomass_lit_sem_L,
                        BAC_TURNOVER_PERC_max, BAC_TURNOVER_PERC_min, GRAZING_EFFECT_hr_max, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max),
               names_to = "SEM_variable", values_to = "SEM") %>% 
  data.frame

biomass_rate_plot$NAME_ORDER <- factor(biomass_rate_plot$NAME, levels = c("Background","Plume", "Quakeplume", "Shrimpocalypse", "Lots 'O Shrimp", "X-18", "Old Man Tree", "Ravelin #2", "Mustard Stand", "Shrimp Hole", "Hot Chimlet #1", "Shrimp Gulley", "South of Hot Chimlet", "South of LungSnack", "Arrow Loop", "Bartizan", "Ravelin #1"))

biomass_rate_plot$VARIABLE_ORDER <- factor(biomass_rate_plot$Variable, 
                                           levels = c("PROK_ml", "EUK_ml",
                                                      "Prok_biomass_L", "Euk_biomass_Hu_L", "Euk_biomass_lit_L",
                                                      "GRAZING_EFFECT_hr", "BAC_TURNOVER_PERC",
                                                      "GRAZE_RATE", "RATE_ugCbac_pergrazer_perday",
                                                      "SPEC_INGESTION_RATE", "SPEC_INGESTION_RATE_lit"),
                                           labels = c("Prokaryote~cells~mL^{-1}", "Eukaryote~cells~mL^{-1}",
                                                      "Prokaryote~µg~C~L^{-1}", "Measured~eukaryote~µg~C~L^{-1}", "Literature-based~eukaryote~µg~C~L^{-1}",
                                                      "Cells~mL^{-1}~hr^{-1}", "Prokaryote~turnover~'%'~d^{-1}",
                                                      "FLP~consumed~min^{-1}", "µg~C~grazer^{-1}~day^{-1}",
                                                      "Cell~carbon~%~day^{-1}", "Cell~carbon~%~day^{-1}~(lit)"))

# View(biomass_rate_plot)

Function to generate consistent plots for Mid-Cayman Rise cell concentration and biomass data.

conc_rate_plot_mcr <- function(df, var, sem){
  df %>% 
    filter(Variable == var) %>%
    filter(SEM_variable == sem) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
    geom_errorbar(aes(ymax = (Value + SEM), ymin = (Value - SEM)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    scale_y_log10() +
    # scale_y_log10(labels = function(x) format(x, scientific = TRUE)) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 14),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")
}
# svg("cells-per-ml-yaxies.svg", h = 8, w = 13)
plot_grid(conc_rate_plot_mcr(biomass_rate_plot, "PROK_ml", "PROK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "EUK_ml", "EUK_sem"),
          conc_rate_plot_mcr(biomass_rate_plot, "Prok_biomass_L", "Prok_biomass_sem_L"),
          conc_rate_plot_mcr(biomass_rate_plot, "Euk_biomass_Hu_L", "Euk_biomass_Hu_sem_L"))

# dev.off()
# svg("biomass-literature-eukC.svg", h=5, w=7)
# conc_rate_plot_mcr(biomass_rate_plot, "Euk_biomass_lit", "Euk_biomass_lit_sem")
# dev.off()
# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
grazing_rate <- biomass_rate_plot %>% 
    filter(Variable == "GRAZING_EFFECT_hr") %>%
    filter(SEM_variable == "GRAZING_EFFECT_hr_min" | SEM_variable == "GRAZING_EFFECT_hr_max") %>% 
    pivot_wider(names_from = SEM_variable, values_from = SEM) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, shape = EXP, fill = SITE)) +
    geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(shape = EXP, fill = SITE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")

# head(biomass_rate_plot)
# svg("plot-grazing-withscale.svg", h = 7, w = 6)
plot_grid(grazing_rate + scale_y_continuous(limits = c(0,2000)),
          grazing_rate,
          ncol = 1, rel_heights = c(0.6,1))
## Warning: Removed 4 rows containing missing values (geom_point).

# dev.off()

# svg("plot-grazing-log.svg", h = 4, w = 6)
# grazing_rate + scale_y_log10(limits = c(1e-1,1e5), breaks=10^(0:7))
grazing_rate + scale_y_continuous(trans=scales::pseudo_log_trans(base = 10),breaks=10^(0:7))

# dev.off()
# scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
prok_turnover <- biomass_rate_plot %>% 
    filter(Variable == "BAC_TURNOVER_PERC") %>%
    filter(SEM_variable == "BAC_TURNOVER_PERC_min" | SEM_variable == "BAC_TURNOVER_PERC_max") %>% 
    pivot_wider(names_from = SEM_variable, values_from = SEM) %>% 
    ggplot(aes(y = Value, x = NAME_ORDER, fill = SITE)) +
    geom_bar(stat = "identity", aes(fill = SITE),
               color = "black", width = 2, position = position_dodge2(preserve = "single")) +
    geom_errorbar(aes(ymax = (BAC_TURNOVER_PERC_max), ymin = (BAC_TURNOVER_PERC_min)), 
                  width = 0.3, position = position_dodge2(preserve = "single")) +
    scale_fill_manual(values = c("#de2d26", "#1c9099")) +
    facet_wrap(VARIABLE_ORDER ~ ., scales = "free", 
               strip.position = c("left"), labeller = label_parsed) +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = "")
# svg("prok-turnover.svg", h = 8, w = 6)
plot_grid(prok_turnover,
          prok_turnover + scale_y_continuous(limits = c(0,100)),
          ncol = 1)
## Warning: Removed 3 rows containing missing values (geom_bar).

# dev.off()

7.6 Compare MCR data with Gorda Ridge - think about environmental parameters

Import Gorda Ridge data

gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/Grazing-calc-wCarbon-results.txt")
# env_gr <- read.delim("../../GordaRidgeCruise_2019/protist-gordaridge-2021/Grazing-at-GordaRidge-SKH-2021/data-input/GR-environ-SAMPLE.txt")

temps <- read.delim("/Users/sarahhu/Desktop/Projects/GordaRidgeCruise_2019/protist-gordaridge-2019-analysis/temperature-allvents.txt")
mcr_graze <- read.delim("table-wcalc.txt")
gr
##              SAMPLE SampleOrigin Bottle    Vent.name  SAMPLE_ORDER       BOTTLE
## 1   BW-Near vent BW           BW    Exp Near vent BW  Near vent BW Experimental
## 2   Vent-Mt Edwards         Vent    Exp   Mt Edwards   Mt. Edwards Experimental
## 3  Vent-Venti latte         Vent    Exp  Venti latte   Venti latte Experimental
## 4   Vent-Candelabra         Vent    Exp   Candelabra    Candelabra Experimental
## 5 Vent-SirVentsalot         Vent    Exp SirVentsalot Sir Ventsalot Experimental
##   cellmL_T0 cellmL_Tf    sem_T0   sem_Tf Hrs_Tf Sample.location  prok_avg
## 1  69818.95  57808.35 6594.5107  472.550     35              BW  51959.11
## 2  60381.13  40639.10 4786.6367 4095.400     24            Vent  51439.52
## 3  45154.57  32448.30  670.3404 5670.600     19            Vent 111192.50
## 4  25622.60  14701.50 3964.4107 3079.078     26            Vent  55076.66
## 5  27985.33  10606.10 5055.5143 4108.875     18            Vent  52998.29
##   MORTALITY MORTALITY_min MORTALITY_max         G     G_min    G_max
## 1 0.1294438    0.06703935     0.1857493  8938.263  4839.402 12329.66
## 2 0.3959460    0.41957541     0.3762200 16818.511 17626.994 16128.80
## 3 0.4174021    0.64113500     0.2325695 31289.007 44259.125 18698.35
## 4 0.5127925    0.57456692     0.4700572 23475.280 25520.903 21977.85
## 5 1.2936683    1.68141744     1.0785052 32912.588 37981.088 29395.13
##   GrazingRate_hr GrazingRate_hr_min GrazingRate_hr_max Prok_turnover
## 1       255.3789           138.2686           352.2759      17.20249
## 2       700.7713           734.4581           672.0331      32.69570
## 3      1646.7899          2329.4276           984.1237      28.13949
## 4       902.8954           981.5732           845.3019      42.62292
## 5      1828.4771          2110.0604          1633.0627      62.10122
##   Prok_turnover_min Prok_turnover_max     N_avg    F_avg         q    G_II_a
## 1          9.313866          23.72954  51959.11 63813.65 0.1882137  9779.414
## 2         34.267415          31.35487  51439.52 50510.12 0.3908531 20105.294
## 3         39.804056          16.81620 111192.50 38801.43 0.3274690 36412.097
## 4         46.337057          39.90411  55076.66 20162.05 0.5416662 29833.162
## 5         71.664736          55.46429  52998.29 19295.72 0.9006783 47734.414
##      G_II_b GrazingRate_hr_II fgC_cell_morono fgC_cell_mcnic
## 1  9779.414          279.4118              86            173
## 2 20105.294          837.7206              86            173
## 3 36412.097         1916.4262              86            173
## 4 29833.162         1147.4293              86            173
## 5 47734.414         2651.9119              86            173
##   cells_consumed_perday fgC_ml_perday_morono fgC_ml_perday_mcnic
## 1              6129.094             527102.1             1060333
## 2             16818.511            1446391.9             2909602
## 3             39522.957            3398974.3             6837471
## 4             21669.489            1863576.1             3748822
## 5             43883.450            3773976.7             7591837
##   ugC_L_perday_morono ugC_L_perday_mcnic lower_mcnichol_morono
## 1           0.5271021           1.060333              3.046833
## 2           1.4463919           2.909602              8.360647
## 3           3.3989743           6.837471             19.647250
## 4           1.8635761           3.748822             10.772116
## 5           3.7739767           7.591837             21.814894
##   upper_mcnichol_morono lower_mcnichol_mcnic upper_mcnichol_mcnic
## 1             0.1640019             6.129094            0.3299108
## 2             0.4500286            16.818511            0.9052901
## 3             1.0575527            39.522957            2.1274025
## 4             0.5798308            21.669489            1.1664037
## 5             1.1742305            43.883450            2.3621148
all_vents <- mcr_graze %>% 
  select(SITE, NAME, SAMPLE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, GRAZING_EFFECT_hr_min, GRAZING_EFFECT_hr_max, BAC_TURNOVER_PERC, ugC_L_perday) %>% 
  rbind(gr %>%
          add_column(SITE = "Gorda Ridge") %>% 
          add_column(EUK_ml = NA) %>% 
          separate(SAMPLE, c("SAMPLE", "NAME"), sep = "-") %>% 
          select(SITE, NAME, SAMPLE, EXP = Bottle, PROK_ml = prok_avg, EUK_ml, GRAZING_EFFECT_hr = GrazingRate_hr, GRAZING_EFFECT_hr_min = GrazingRate_hr_min, GRAZING_EFFECT_hr_max = GrazingRate_hr_max, BAC_TURNOVER_PERC = Prok_turnover, ugC_L_perday = ugC_L_perday_morono)) %>% 
  left_join(temps) %>% 
  mutate(SAMPLE_TYPE = case_when(
    grepl("BSW", NAME) ~ "Background",
    grepl("Near vent BW", NAME) ~ "Background",
    grepl("Background", NAME) ~ "Background", 
    grepl("Plume", NAME) ~ "Background",
    TRUE ~ "Vent"
  ))
## Joining, by = c("NAME", "SAMPLE")
# View(all_vents)
all_vents
##           SITE           NAME                 SAMPLE EXP   PROK_ml    EUK_ml
## 1     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 2     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 3     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 4     Von Damm          Plume               VD-Plume Bag  16478.31 157.77468
## 5     Von Damm           X-18                 VD-X18 Bag 111429.78 314.87222
## 6     Von Damm   Old Man Tree                 VD-OMT IGT  71147.40 472.30833
## 7     Von Damm     Ravelin #2                VD-Rav2 Bag  71147.40 409.33389
## 8     Von Damm     Ravelin #2                VD-Rav2 IGT  71147.40 620.99799
## 9     Von Damm     Ravelin #2                VD-Rav2 IGT  71147.40 620.99799
## 10    Von Damm  Mustard Stand        VD-MustardStand Bag  56677.00 259.76958
## 11    Von Damm    Shrimp Hole          VD-ShrimpHole Bag  41982.96 385.71847
## 12     Piccard          Plume          Piccard-Plume Bag  51429.13  79.30115
## 13     Piccard Shrimpocalypse Piccard-Shrimpocalypse Bag 238585.68 454.81543
## 14     Piccard Shrimpocalypse Piccard-Shrimpocalypse IGT 238585.68 454.81543
## 15     Piccard Lots 'O Shrimp    Piccard-LotsOShrimp Bag  53878.14 230.90630
## 16     Piccard Lots 'O Shrimp    Piccard-LotsOShrimp Bag  53878.14 230.90630
## 17 Gorda Ridge   Near vent BW                     BW Exp  51959.11        NA
## 18 Gorda Ridge     Mt Edwards                   Vent Exp  51439.52        NA
## 19 Gorda Ridge    Venti latte                   Vent Exp 111192.50        NA
## 20 Gorda Ridge     Candelabra                   Vent Exp  55076.66        NA
## 21 Gorda Ridge   SirVentsalot                   Vent Exp  52998.29        NA
##    GRAZING_EFFECT_hr GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max
## 1          127.06388              96.81058             157.31719
## 2          127.06388              96.81058             157.31719
## 3          127.06388              96.81058             157.31719
## 4           24.02685              13.80868              34.24501
## 5         1166.28778             777.52519            1555.05037
## 6         5465.73080            4048.68948            6882.77212
## 7          116.85656              95.88230             137.83081
## 8         7624.36451            6112.04639            9136.68264
## 9          793.88963             636.41897             951.36028
## 10           0.00000               0.00000               0.00000
## 11           0.00000               0.00000               0.00000
## 12          44.25930              34.87229              53.64630
## 13        6006.74315                    NA                    NA
## 14       20427.17803           17284.53525           23569.82080
## 15           0.00000                    NA                    NA
## 16           0.00000                    NA                    NA
## 17         255.37893             138.26863             352.27592
## 18         700.77129             734.45809             672.03314
## 19        1646.78986            2329.42764             984.12374
## 20         902.89538             981.57320             845.30189
## 21        1828.47708            2110.06045            1633.06274
##    BAC_TURNOVER_PERC ugC_L_perday Highest.Temp  pH   Mg  H2S    H2  CH4
## 1           8.048465   0.26225986        4.181  NA   NA   NA    NA   NA
## 2           8.048465   0.26225986        4.100  NA   NA   NA    NA   NA
## 3           8.048465   0.26225986        4.000  NA   NA   NA    NA   NA
## 4           3.499414   0.04959142        4.230  NA   NA   NA    NA   NA
## 5          25.119772   2.40721798       62.000  NA   NA   NA    NA   NA
## 6         184.374334  11.28126837      121.000  NA   NA   NA    NA   NA
## 7           3.941897   0.24119193      112.000  NA   NA   NA    NA   NA
## 8         257.191065  15.73668836      112.000  NA   NA   NA    NA   NA
## 9          26.780110   1.63858819      112.000  NA   NA   NA    NA   NA
## 10          0.000000   0.00000000      108.000  NA   NA   NA    NA   NA
## 11          0.000000   0.00000000       22.000  NA   NA   NA    NA   NA
## 12          2.065411   0.09135119        4.722  NA   NA   NA    NA   NA
## 13         60.423507  12.39791787       90.000  NA   NA   NA    NA   NA
## 14        205.482690  42.16169545       90.000  NA   NA   NA    NA   NA
## 15          0.000000   0.00000000       36.000  NA   NA   NA    NA   NA
## 16          0.000000   0.00000000       35.000  NA   NA   NA    NA   NA
## 17         17.202493   0.52710212           NA  NA   NA   NA    NA   NA
## 18         32.695699   1.44639193       30.000 5.8 42.8 1.01 127.0 10.1
## 19         28.139494   3.39897427       23.000 5.5 50.4   NA    NA  0.9
## 20         42.622919   1.86357606       68.000 5.8 45.8   NA  21.9 23.7
## 21         62.101220   3.77397670       53.000  NA 50.8   NA    NA   NA
##    SAMPLE_TYPE
## 1   Background
## 2   Background
## 3   Background
## 4   Background
## 5         Vent
## 6         Vent
## 7         Vent
## 8         Vent
## 9         Vent
## 10        Vent
## 11        Vent
## 12  Background
## 13        Vent
## 14        Vent
## 15        Vent
## 16        Vent
## 17  Background
## 18        Vent
## 19        Vent
## 20        Vent
## 21        Vent

Plot grazing rates

allrates <- all_vents %>% 
  select(SITE, NAME, SAMPLE, EXP, SAMPLE_TYPE, starts_with("GRAZING_EFFECT_")) %>% 
  distinct() %>% 
    ggplot(aes(y = GRAZING_EFFECT_hr, x = NAME, fill = SITE, shape = SAMPLE_TYPE)) +
    geom_errorbar(aes(ymax = (GRAZING_EFFECT_hr_max), ymin = (GRAZING_EFFECT_hr_min)), 
                  width = 0.2, position = position_dodge(width = 0.4)) +
    geom_point(stat = "identity", aes(fill = SITE, shape = SAMPLE_TYPE),
               color = "black", size = 3, position = position_dodge(width = 0.4)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#de2d26", "#1c9099", "#addd8e")) +
    facet_grid(. ~ SAMPLE_TYPE, scales = "free", space = "free") +
    theme_minimal() +
    theme(panel.grid.major = element_line(), panel.grid.minor = element_blank(),
          panel.background = element_blank(), 
          axis.line = element_line(colour = "black"), 
          axis.text.x = element_text(color="black", size = 11, 
                                     angle = 45, hjust = 1, vjust = 1), 
          axis.text.y = element_text(color="black", size = 11),
          axis.title =element_text(color="black", size = 11),
          axis.ticks = element_line(),
          legend.title = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(color="black", size = 11),
          strip.text.x = element_blank())+
    guides(fill = guide_legend(override.aes = list(shape = c(21))),
           shape = guide_legend(override.aes = list(fill = "black"))) +
    labs(x = "", y = bquote("cells"~ml^-1~hr^-1))
# svg("compare-all-rates.svg", h = 10, w = 7)
# plot_grid(allrates,
          # allrates + scale_y_continuous(limits = c(0,1000)),
          # allrates + scale_y_log10(), ncol = 1)
# dev.off()
# svg("compare-all-rates-color.svg", h = 4, w = 9)
allrates + scale_y_log10()
## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

# dev.off()

7.7 Run regression analysis with MCR and GR data

metadata <- read.delim("input-data/flp-exp-metadata-compiled.txt")
# head(metadata)
# temps
all_vents
##           SITE           NAME                 SAMPLE EXP   PROK_ml    EUK_ml
## 1     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 2     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 3     Von Damm     Background                 VD-BSW Bag  37889.62  91.83773
## 4     Von Damm          Plume               VD-Plume Bag  16478.31 157.77468
## 5     Von Damm           X-18                 VD-X18 Bag 111429.78 314.87222
## 6     Von Damm   Old Man Tree                 VD-OMT IGT  71147.40 472.30833
## 7     Von Damm     Ravelin #2                VD-Rav2 Bag  71147.40 409.33389
## 8     Von Damm     Ravelin #2                VD-Rav2 IGT  71147.40 620.99799
## 9     Von Damm     Ravelin #2                VD-Rav2 IGT  71147.40 620.99799
## 10    Von Damm  Mustard Stand        VD-MustardStand Bag  56677.00 259.76958
## 11    Von Damm    Shrimp Hole          VD-ShrimpHole Bag  41982.96 385.71847
## 12     Piccard          Plume          Piccard-Plume Bag  51429.13  79.30115
## 13     Piccard Shrimpocalypse Piccard-Shrimpocalypse Bag 238585.68 454.81543
## 14     Piccard Shrimpocalypse Piccard-Shrimpocalypse IGT 238585.68 454.81543
## 15     Piccard Lots 'O Shrimp    Piccard-LotsOShrimp Bag  53878.14 230.90630
## 16     Piccard Lots 'O Shrimp    Piccard-LotsOShrimp Bag  53878.14 230.90630
## 17 Gorda Ridge   Near vent BW                     BW Exp  51959.11        NA
## 18 Gorda Ridge     Mt Edwards                   Vent Exp  51439.52        NA
## 19 Gorda Ridge    Venti latte                   Vent Exp 111192.50        NA
## 20 Gorda Ridge     Candelabra                   Vent Exp  55076.66        NA
## 21 Gorda Ridge   SirVentsalot                   Vent Exp  52998.29        NA
##    GRAZING_EFFECT_hr GRAZING_EFFECT_hr_min GRAZING_EFFECT_hr_max
## 1          127.06388              96.81058             157.31719
## 2          127.06388              96.81058             157.31719
## 3          127.06388              96.81058             157.31719
## 4           24.02685              13.80868              34.24501
## 5         1166.28778             777.52519            1555.05037
## 6         5465.73080            4048.68948            6882.77212
## 7          116.85656              95.88230             137.83081
## 8         7624.36451            6112.04639            9136.68264
## 9          793.88963             636.41897             951.36028
## 10           0.00000               0.00000               0.00000
## 11           0.00000               0.00000               0.00000
## 12          44.25930              34.87229              53.64630
## 13        6006.74315                    NA                    NA
## 14       20427.17803           17284.53525           23569.82080
## 15           0.00000                    NA                    NA
## 16           0.00000                    NA                    NA
## 17         255.37893             138.26863             352.27592
## 18         700.77129             734.45809             672.03314
## 19        1646.78986            2329.42764             984.12374
## 20         902.89538             981.57320             845.30189
## 21        1828.47708            2110.06045            1633.06274
##    BAC_TURNOVER_PERC ugC_L_perday Highest.Temp  pH   Mg  H2S    H2  CH4
## 1           8.048465   0.26225986        4.181  NA   NA   NA    NA   NA
## 2           8.048465   0.26225986        4.100  NA   NA   NA    NA   NA
## 3           8.048465   0.26225986        4.000  NA   NA   NA    NA   NA
## 4           3.499414   0.04959142        4.230  NA   NA   NA    NA   NA
## 5          25.119772   2.40721798       62.000  NA   NA   NA    NA   NA
## 6         184.374334  11.28126837      121.000  NA   NA   NA    NA   NA
## 7           3.941897   0.24119193      112.000  NA   NA   NA    NA   NA
## 8         257.191065  15.73668836      112.000  NA   NA   NA    NA   NA
## 9          26.780110   1.63858819      112.000  NA   NA   NA    NA   NA
## 10          0.000000   0.00000000      108.000  NA   NA   NA    NA   NA
## 11          0.000000   0.00000000       22.000  NA   NA   NA    NA   NA
## 12          2.065411   0.09135119        4.722  NA   NA   NA    NA   NA
## 13         60.423507  12.39791787       90.000  NA   NA   NA    NA   NA
## 14        205.482690  42.16169545       90.000  NA   NA   NA    NA   NA
## 15          0.000000   0.00000000       36.000  NA   NA   NA    NA   NA
## 16          0.000000   0.00000000       35.000  NA   NA   NA    NA   NA
## 17         17.202493   0.52710212           NA  NA   NA   NA    NA   NA
## 18         32.695699   1.44639193       30.000 5.8 42.8 1.01 127.0 10.1
## 19         28.139494   3.39897427       23.000 5.5 50.4   NA    NA  0.9
## 20         42.622919   1.86357606       68.000 5.8 45.8   NA  21.9 23.7
## 21         62.101220   3.77397670       53.000  NA 50.8   NA    NA   NA
##    SAMPLE_TYPE
## 1   Background
## 2   Background
## 3   Background
## 4   Background
## 5         Vent
## 6         Vent
## 7         Vent
## 8         Vent
## 9         Vent
## 10        Vent
## 11        Vent
## 12  Background
## 13        Vent
## 14        Vent
## 15        Vent
## 16        Vent
## 17  Background
## 18        Vent
## 19        Vent
## 20        Vent
## 21        Vent

Linear regression with both Gorda Ridge and MCR data

library(broom)
# ?pivot_longer
regression_input <- all_vents %>% 
  filter(!is.na(Highest.Temp)) %>% 
    select(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, PROK_ml, EUK_ml, GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday, TEMP = Highest.Temp) %>% 
  pivot_longer(cols = c(GRAZING_EFFECT_hr, BAC_TURNOVER_PERC, ugC_L_perday), names_to = "RATE", values_to = "RATE_VALUE") %>% 
  pivot_longer(cols = c(PROK_ml, EUK_ml, TEMP), names_to = "PARAMS", values_to = "PARAMS_VALUE") %>% 
  data.frame

regression_tmp <- regression_input %>% 
  # Set up the linear regression
  group_by(RATE, PARAMS) %>% 
  nest(-RATE, -PARAMS) %>% 
  mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)), 
    tidied = map(lm_fit, tidy)) %>% 
  unnest(tidied) %>% 
  select(RATE, PARAMS, 
    term, estimate) %>% 
  pivot_wider(names_from = term, values_from = estimate) %>% 
  select(everything(), SLOPE = PARAMS_VALUE) %>%
  data.frame
## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?
regression_results <- regression_input %>% 
  group_by(RATE, PARAMS) %>% 
  nest(-RATE, -PARAMS) %>% 
  mutate(lm_fit = map(data, ~lm(RATE_VALUE ~ PARAMS_VALUE, data = .)),
    glanced = map(lm_fit, glance)) %>% 
  unnest(glanced) %>% 
  select(RATE, PARAMS, r.squared, adj.r.squared) %>% 
  right_join(regression_tmp) %>%
  right_join(regression_input) %>%
  data.frame
## Warning: All elements of `...` must be named.
## Did you want `data = c(SITE, NAME, SAMPLE, SAMPLE_TYPE, EXP, RATE_VALUE, PARAMS_VALUE)`?
## Joining, by = c("RATE", "PARAMS")
## Joining, by = c("RATE", "PARAMS")
# head(regression_results)

Plot regression results

regression_results %>% 
  # filter(RATE == "GRAZING_EFFECT_hr") %>%
  # filter(PARAMS == "TEMP") %>% 
  ggplot(aes(x = PARAMS_VALUE, y = RATE_VALUE, shape = SAMPLE_TYPE, fill = SITE)) +
    geom_abline(aes(slope = SLOPE, intercept = `X.Intercept.`), color = "black", linetype = "dashed", size = 0.5) +
    geom_point(color = "black", aes(shape = SAMPLE_TYPE, fill = SITE)) +
    scale_shape_manual(values = c(21, 24)) +
    scale_fill_manual(values = c("#476AA7","#7299CE", "#A2937A")) +
    facet_wrap(PARAMS ~ RATE + round(r.squared, 3), scales = "free", ncol = 5,
             strip.position = "bottom", labeller = label_parsed) +
  theme_bw() +
  theme(
    strip.background = element_blank(),
    strip.placement = "outside",
    strip.text = element_text(color = "black", size = 10),
    axis.title = element_text(color = "black", size = 10),
    legend.title = element_blank()) # +
## Warning: Removed 12 rows containing missing values (geom_point).

  # labs(y = bquote("Cells "~mL^-1 ~hr^-1), x = "Temperature (C)")
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS/LAPACK: /Users/sarahhu/anaconda3/envs/r_3.6.0/lib/R/lib/libRblas.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] broom_0.7.5     cowplot_1.0.0   forcats_0.5.0   stringr_1.4.0  
##  [5] dplyr_1.0.4     purrr_0.3.4     readr_1.3.1     tidyr_1.1.3    
##  [9] tibble_3.1.0    ggplot2_3.3.1   tidyverse_1.3.0
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.20          haven_2.3.1        colorspace_2.0-0  
##  [5] vctrs_0.3.6        generics_0.1.0     htmltools_0.5.1.1  yaml_2.2.1        
##  [9] utf8_1.1.4         blob_1.2.1         rlang_0.4.10       pillar_1.5.0      
## [13] withr_2.4.1        glue_1.4.2         DBI_1.1.0          RColorBrewer_1.1-2
## [17] dbplyr_1.4.4       modelr_0.1.8       readxl_1.3.1       lifecycle_1.0.0   
## [21] munsell_0.5.0      gtable_0.3.0       cellranger_1.1.0   rvest_0.3.5       
## [25] evaluate_0.14      labeling_0.4.2     knitr_1.31         fansi_0.4.2       
## [29] highr_0.8          Rcpp_1.0.5         scales_1.1.1       backports_1.2.1   
## [33] debugme_1.1.0      jsonlite_1.6.1     farver_2.1.0       fs_1.4.1          
## [37] hms_0.5.3          digest_0.6.27      stringi_1.4.6      grid_3.6.1        
## [41] cli_2.3.1          tools_3.6.1        magrittr_2.0.1     crayon_1.4.1      
## [45] pkgconfig_2.0.3    ellipsis_0.3.1     xml2_1.3.2         reprex_0.3.0      
## [49] lubridate_1.7.8    assertthat_0.2.1   rmarkdown_2.6      httr_1.4.2        
## [53] rstudioapi_0.11    R6_2.5.0           compiler_3.6.1